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Fourth - order gravity tiieories have received much interest in recent years thanks to their ability 
to provide an accelerated cosmic expansion in a matter only universe. In these theories, the La- 
grangian density of the gravitational field has the form R + f{R), and the explicit choice of the 
arbitrary function f{R) must meet the local tests of gravity and the constraints from the primordial 
abundance of the light elements. Two popular classes of f{R) models, which are expected to fulfill 
all the above requirements, have recently been proposed. Ifowever, neither of these models has ever 
been quantitatively tested against the available astrophysical data. Here, by combining Type la 
Supernovae and Gamma Ray Bursts, we investigate the ability of these models to reproduce the 
observed Hubble diagram over the redshift range (0, 7). We find that both models fit very well this 
dataset with the present-day values of the matter density and deceleration parameters which agree 
with previous estimates. However, the strong degeneracy among the f{R) parameters prevents us 
from putting strong constraints on the values of these parameters; nevertheless, we can identify the 
regions of the parameter space that should, in principle, be carefully explored with future data and 
dynamical probes in order to discriminate among f{R) theories and standard dark energy models. 

PACS numbers: 98.80.-k, 98.80.Es, 95.36.+d, 95.36.+X 



I. INTRODUCTION 

It is no-w -widely accepted that the universe is presently 
undergoing a period of accelerated expansion. The Hub- 
ble diagram of Type la Supernovae (SNela) H i, i, i, 
H, @] J the anisotropy spectrum of the cosmic micro-wave 
background radiation (CMBR) 0, H, Q and the cluster- 
ing properties probing the large-scale structure of the 
universe IflJ are concordant pieces of evidence in favour 
of this cosmic speed up. In the concordance ACDM cos- 
mological model, this -wide dataset is excellently repro- 
duced [T]| by assuming a spatially flat universe domi- 
nated by Cold Dark Matter (CDM) and a cosmological 
constant A [l^ . Ho-wever, this scenario has t-wo serious 
dra-wbacks: (1) if interpreted as vacuum energy, the A 
value is 120 orders of magnitude smaller than -what is 
expected from quantum field theory; (2) the coincidence 
and fine - tuning problems do not seem to have a natural 
explanation. This circumstance has motivated the search 
for alternative explanations. They mostly rely on scalar 
fields -with a suitable potential -which provides a vary- 
ing A term or other dark energy fluids 13] -with exotic 
properties. 

An alternative route towards solving the problem of 
the accelerated cosmic expansion may be explored by 
looking at the cosmological constant as an additive con- 
stant term in the Einstein- Hilbert gravity Lagrangian: 
from the point of vie-w of the field equations, one is not 
adding any ne-w source term to the energy-momentum 
tensor, but rather it is the geometrical sector (the left- 
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hand side of the equations) that has been modified. It 
is therefore worth wondering whether the accelerated ex- 
pansion can indeed be obtained by generalizing this ap- 
proach, i.e. by adding other non- linear terms to the 
action. The search for a generalization of Einstein Gen- 
eral Relativity (GR) actually dates back to just few years 
after the seminal Einstein papers (see, e.g., [l3| for a 
historical review), but the relevance of these corrective 
terms were considered to be restricted only to the strong 
gravity regimes. Moreover, effective field theory consid- 
erations motivated the introduction of small couplings 
which suppress these non - linear terms in all the other 
curvature regimes. As a consequence, corrections to GR 
were considered only close to the Planck scale, i.e. in 
the very early - universe llSlI . and in the attempt to avoid 
black hole singularities [iq . 

The qualitative similarities between the present - day 
acceleration and the early universe inflationary expansion 
renewed the interest in what are now collectively referred 
to as f{R) theories. Indeed, it became clear soon that 
an accelerated expansion in a matter only universe can 
be achieved by suitably choosing the functional expres- 
sion for f{R) leading to a non-linear gravity Lagr ang ian, 
C cx R. + f{R). The impressive amount of papers [iTlflsj 
investigating different aspects of f{R) theories helps un- 
derstanding that the success in explaining the cosmic ac- 
celeration must be balanced by the constraints provided 
by the local gravity tests (see, e.g., [l^ for exhaustive 
reviews discussing these issues). Indeed, many models 
providing an accelerated expansion were later rejected 
because the non-linear terms do not turn off their effect 
on the Solar System scale and thus lead to unacceptable 
contrast with the success of GR in this regime. 

Recently, two models carefully designed to evade the 
local gravity tests but still providing an accelerated cos- 
mic expansion has been proposed [20. [28|. However, nei- 
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ther of these models has yet been quantitatively tested 
against the available astrophysical data. Here, in order to 
to further support (or rule out) f{R) theories as alterna- 
tive candidates to the dark energy hypothesis, we search 
for the appropriate set of parameters of these models to 
verify that the acceleration, that they qualitatively pre- 
dict, actually quantitevely agree with the data. 

The scheme of the paper is as follows. In § II, we will 
briefly remind the basics of f{R) theories writing down 
the equations needed to determine the background evo- 
lution of the universe. § III discusses the theoretical con- 
straints driving the choice of a functional expression for 
f{R) thus motivating the adoption of the two popular 
models to be tested. The data used and the statistical 
method used to constrain the model parameters are pre- 
sented in § IV, while the results of the analysis are com- 
mented upon in § V. Conclusions and future perspectives 
are finally given in § VI. 



II. f{R) COSMOLOGY 

The idea of modifying the gravity sector of the Ein- 
steinian General Relativity (GR) dates back to the 
Starobinsky [l5l | attempts to get an inflationary expan- 
sion of the early universe without the need of any scalar 
field. This can be accomplished by generalizing the Ein- 
stein - Hilbert action as : 
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where R is the scalar curvature, = SttG (we use units 
where the light speed c = 1), Cm is the standard matter 
Lagrangian, and f{R) is an analytical function expressing 
the deviation from the Einstein GR. For f{R) = — 2A, we 
get back the concordance ACDM model, while nontrivial 
dynamics is obtained for other choices. 

Under the metric approach to f{R) gravity, the field 
equations are obtained by varying the action ([T]) with 
respect to the metric only. We obtain : 

Ga(i + fRRafl - (^f - n/fl,^ gap - = l^'^Tafj 

(2) 

with Gap = Rap — {R/2)gap and Tap the usual Einstein 
and source stress - energy tensor. Hereafter the subscript 
R will denote differentiation with respect to R. We also 
assume a spatially flat Robertson - Walker metric, with 
scale factor a; in this case, the scalar curvature reads 



R = UH^ + 6HH' 



(3) 



where H = d/a is the Hubble parameter and the dot in- 
dicates the time derivative; by using the more convenient 
variable ?7 = Ina, we denote the derivative with respect 
to 1] with a prime. With these definitions, Eqs.([2]) lead 



to the modified Friedmann equation 

^ (HH' + H^) + 1/ + H'fRRR' - ypM , (4) 

where we have assumed that dust matter with enery den- 
sity pm and pressure pm = is the only fluid filling the 
universe. To solve Eqs.Q and ([4]), we follow [20| and 
introduce the dimensionless variables : 



VH = — 
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with 
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a convenient curvature scale depending on the present 
day values of the matter density parameter Qm and the 
Hubble constant h = iJo/(100 km/s/Mpc). The back- 
ground evolution of the universe is determined by the 
following set of coupled ordinary differential equations : 
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(7) 
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It is worth noting that the above system of equations 
of first order in (y_f/,y/j) is equivalent to a single second 
order equation in yfj. If we look at the definition of yn, 
we see that the full system is equivalent to a single fourth- 
order non-linear differential equation for the scale factor 
a(t): this property explains why /(i?) theories are also 
commonly referred to as fourth-order gravity theories. 

In order to integrate the system, we need to set the 
boundary conditions. In |20i], the authors set them at red- 
shifts z — > oo, by requiring that f{R) tends to a constant 
at this epoch and by considering the detailed balance 
of the perturbative corrections to i? = k^Pm- However, 
here we are interested in fitting the model to data that 
probe the range z < 10; therefore, it is more interesting 
to set the present day values of {vhtUr) and integrate 
the equations back in time up to a ~ 0.001 (77 ~ By 
setting, as usual, ao = 1 and remembering that 

R^ = ml[l-qo) (9) 

with q ~ —aa/b? the deceleration parameter, we get: 

VH[0) = Hl/m'-l , yR{Q)^&{Hl/m'){l-qa)-i . 

(10) 

It is worth noting that, because of Eq.®, the initial con- 
ditions are determined by the values of the three pa- 
rameters (r^M, qo). This is again consistent with f{R) 
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models being fourth-order theories; thus, we require three 
initial conditions to determine the evolution of the scale 
factor for any given f{R). Had we written down a sin- 
gle equation for a{t), we should have set the present-day 
values of the first three time derivatives of a(t) which in- 
clude the jerk parameter j — {<Pa/dt^)H^^ . In contrast, 
the introduction of the (j/HtUr) auxiliary functions and 
of the curvature scale m enables us to replace the quite 
uncertain jerk parameter with the more manageable mat- 
ter density ^Im- 



III. f{R) MODELS 



A key role in fourth - order gravity is obviosuly played 
by the functional expression of f{R). In principle, such 
a choice is fully arbitrary unless one has a theoretical 
motivation leading to a unique expression for f{R). Ac- 
tually, things are not so easy. Indeed, the modification 
of GR introduces deviations not only on the cosmological 
scales, but at all the scales where gravitational phenom- 
ena can be tested. In particular, it has proven to be 
quite difficult for a large class of f{R) theories to evade 



the constraints on the Solar System scale (see, e.g., [21 1 
and refs therein). As Chiba [22] has shown, the main dif- 
ficulty arises from /(i?) models introducing a new scalar 
degree of freedom with the same coupling to matter as 
gravity. As a consequence, it appears a long range fifth 
force that violates the constraints on the PPN parame- 
ters. Although the derivation of the PPN parameters has 
been questioned [l^ , significant deviations from the GR 
metric around the Sun seem to be confirmed because of 
a decoupling of the scalar curvature from the local den- 
sity. As a possible way out, one can invoke a chamaleon 
effect [2J] to reassociate high density with high curvature 
so that the scalar degree of freedom becomes very mas- 
sive and the fifth force escapes any detection. To this 
aim, f{R) should be tailored in such a way to give rise 
to a mass squared term which is large and positive in 
high curvature environments [2^. It is worth stressing 
that this same condition is also required if we want to re- 
cover GR in the early universe [20| and obtain the usual 
matter dominated era. Since i? ^ oo in this limit, we 
expect that f{R) tends to a small constant in order to 
make its effect negligible with respect to the GR term. 
On the other hand, in the late universe, we expect to 
mimic the same evolutionary history of the ACDM be- 
cause this model agrees with the data; therefore f{R) 
should reduce again to a small constant, but it should 
again tend to zero in the limit of a vanishing R to agree 
with the observational fact that A takes a very low value. 
Summarizing, one has to look for a functional expression 



satisfying the following constraints: 



f lim /(i?)=0 

IX — >U 



lim f{R) — const 

i?— i-OO 



, (11) 



fR{R)\ = df{R)/dR\Ryy^2 > 

. fRR{R)\R»m- = d^f{R)/dR^\Ryy^2 > 

where is a typical curvature scale and the last condi- 
tion [13] ensures that, in the limit R ^ m^, the solution 
is stable at high curvature. 

Among the possible choices left out by the above con- 
ditions, we will consider here two classes of f{R) models. 
For the first one, we follow [20] and set: 



fiR) 



ci(i?/m2)" 
TTc^iR/m^ 



(12) 



with m given by ([6|), and (n, ci,C2) are positive dimen- 
sionless constants. We will refer to this choice as the 
Hu & Sawicki (hereafter, HS) model by the name of the 
authors who first suggested this expression. It is easy to 
check that all the constraints (fTTj) are easily passed by the 
HS model. In particular, we note that, since /(O) = 0, 
there is no actual cosmological constant in the model, 
but 



/R~>0 



lim f{R) ~ -^m^ + I — 

" " " C2 " 



Cl 



m' 
R 



so that, when ci/c^ — > at fixed Ci/c2, we recover an ef- 
fective cosmological constant in high curvature (m^ /R 
0) environments. 

Another possibility to satisfy all the constraints (|lip is 
offered by the Starobinsky proposal [28|] : 



f{R) = XR. 



R^ 
Rl 



- 1 



(13) 



with Ri, a scaling curvature parameter and (A, n) two 
positive constants. We will refer to this class of J{R) 
theories as the Starobinsky (St) model. Note that, even 
in this case, /(O) = so that no actual cosmological 
constant is present. Nevertheless, an effective one is re- 
covered in the high curvature regime as can be seen from 
f{R > R^) - -2Aoo with Aoo = \R./2. 

As a general remark, we note that the HS and St mod- 
els are quite similar at both very low and very high red- 
shifts since they are both built up by imposing the same 
constraints on f{R). Moreover, they both aim at mim- 
icking the successful ACDM scenario in the late and early 
universe. Put in other words, the HS and St models both 
reduces to the GR -t- A case in the limits of very high and 
very low curvature. What makes them different is the 
way the two extreme cases are connected, i.e. how the 
universe evolves from the present day A dominated phase 
to the early matter epoch. 
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For completeness, we finally remind the reader that the 
two models we are considering here are not the o nly vi- 
able ones; other possible examples are given in [2^[30]- It 
is, moreover, possible to work out /(i?) models which can 
also provide an inflationary expansion in the very early 
universe [sH- However, all these other cases share many 
similarities with the HS and St models so that we are 
confident that exploring just two classes of fourth-order 
gravity theories should provide us with some general con- 
clusions on their viability. 



IV. CONSTRAINING THE MODELS 

Any model that aims to describe the evolution of the 
universe must be able to reproduce what is indeed ob- 
served. Matching the model with observations is also a 
powerful tool to constrain its parameters and allows to 
estimate some further quantities of interest. 

As mentioned above, the HS and St models are care- 
fully designed to give an effective cosmological constant 
in the late universe. Moreover, both theories are assigned 
by a three parameters function, and one can anticipate 
that the considerable freedom allowed by the degeneracy 
among the parameters makes it easy to find some com- 
binations that lead to quite similar H{z) in the low z 
regime. Fitting to SNela data can only tell us whether 
the models are viable over the redshift range (0, 1.5). We 
expect that this is indeed the case because the f{R) func- 
tions have been tailored to do so. What is not garanteed 
is that the HS and St f{R) models can describe the back- 
ground expansion up to higher redshifts z. To probe this 
regime, we will use the recently derived Hubble diagram 
of Gamma Ray Bursts (GRBs) . 

In the Bayesian approach to model testing, we explore 
the parameter space through the likelihood function : 



£(p) CX CsNelaip) CgB.b{p) X 



X exp 



X exp 



1 /lj' 



1 / husT - h 

2 V (^HST 



(14) 



where p denotes the set of model parameters. Before dis- 
cussing in detail the term related to the SNela and GRB 
data, we concentrate on the two Gaussian priors. The 
former takes into account the constraints on the physical 
matter density lom ~ ^Kjh^ with 



, ,obs _i_ ^ 



0.137 ±0.004 



as inferred from WMAP5 data (9|]. One can wonder 
whether such an estimate may be used as a constraint 
on the /(i?) models since it has been obtained by fitting 
the CMBR spectrum assuming the validity of GR. How- 
ever, what is really needed for this estimate to be model 
independent is not that GR holds along all the evolution- 
ary history, but that the gravity Lagrangian reduces to 



the Einstein- Hilbert one at the last scattering. Since, 
for both HS and St models, f{R)/R ^ for z ~ 1000, 
we can safely use the WMAP5 uim value as a constraint. 

The Gaussian prior on h in Eq. (ll4p stems from the 
results of the HST Project [32| which has estimated the 
Hubble constant Hq using a well calibrated set of local 
distance scale. Averaging over the different methods, the 
survey finally gives : 

hnsT ± (THST = 0.72 ± 0.08 

as a cosmological model independent constraint"'^. 

The main two terms in the likelihood are both re- 
lated to the Hubble diagram, the first one being, in par- 
ticular, connected with SNela. These latter data have 
provided the first piece of evidence for the cosmic speed 
up and are still a sort of ground-zero test that every cos- 
mological model has to pass to be considered acceptable. 
To check this, one relies on the predicted distance mod- 
ulus : 



Hthiz,p) = 25 + Slog 



— {l + z)r{z,p) 

^0 



with r(z) the dimensionless comoving distance : 

The likelihood function is then defined as 

1 



(15) 



(16) 
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SNel 



a(p) 



X exp 



(17) 



where AfsNeia is the total number of SNela used, A/i is a 
■^SNeia - dimensional vector with the values of Hohsizi) — 
fithizi) and CsNeia is the AfsNeia X AfsNeia covariauce 
matrix of the SNela data. Note that, if we neglect the 
correlation induced by systematic errors^, as we do here, 
C SNela is a diagonal matrix and Eq. (jl7|) simplifies to : 



CsNelaip) OC exp [-XsWe/a (P) /2] 



with 



xlNelaiP) = 51 



^■obs{zi) - Hth{zi) 



(18) 



(19) 



While this work was near completion, the SHOES collaboration 
[s^ l has released a more precise estimate as h = 0.742 ± 0.036 in 
agreement with our adopted value. 

Actually, it has been shown |34| that neglecting systematic er- 
rors does not shift the central values, but only weakens the con- 
straints. Although such a result has been obtained using stan- 
dard dark energy models, we are confident that this is also the 
case for our /(-R) theories. 
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with (Ti the error on the observed distance modulus 
Hobsizi) for the i-th object at redshift z^. As input data, 
we use the Union SNela sample assembled in [34| by re- 
analysing with the same pipeline both the recent SNela 
SNLS [J] and ESSENCE [5,] samples and older nearby 
and high redshift (sj datasets. 

Although quite useful in probing the accelerated ex- 
pansion, SNela are limited to z ~ 1.5. As a conse- 
quence, one has to resort to a different distance indi- 
cator to push the Hubble diagram to higher redshift and 
probe the (supposedly) matter dominated era. Thanks to 
the enormous energy release that makes them visibile up 
to z ~ 6.6, GRBs stand out as ideal candidates to this 
scope. The discovery of 2D correlations between their 
properties have opened the way towards making GRBs 
standard candles similarly to SNela [s^. As a result, 
Schaefer [s^ have provided the first GRBs Hubble di- 
agram containing 69 objects with /iobs(-z) estimated by 
averaging over 5 different 2D correlations. We use here 
the updated GRBs Hubble diagram recently presented in 
[3^ based on a model-independent recalibration of the 
same 2D correlations used by Schaefer. 

Since there is no correlation among the errors of dif- 
ferent GRBs, the likelihood function now simply reads : 



Cgrb{p) oc exp [-XGi?i3(p)/2] 



with 



XGRsiP) 



J^GRB 



(20) 



(21) 



where ugrb takes care of the intrinsic scatter inherited 
from the scatter of GRBs around the 2D correlations used 
to derive the individual distance moduli.'^ 




FIG. 1: Best fit HS model superimposed to the data. 



the most plausible model in an Occam's razor sense given 
the data at hand. However, in a Bayesian context, the 
best fit parameters individually do not necessarily have 
to be probable, but rather they must have a high joint 
probability density that might occupy only a small vol- 
ume of the parameter space. This situation can arise if 
the best fit solution does not lie in the bulk of the pos- 
terior probability distribution. Such a situation may of- 
ten occur when the posterior is non-symmetric in a high 
dimensional space so that the volume can dramatically 
increase with the distance from the best fit solution. In 
this case, the best fit solution for each parameter could 
easily lie outside the bulk of the individual posterior dis- 
tribution for Pi obtained by marginalizing over the other 
parameters. This is indeed what happens for our mod- 
els so that we have preferred to remind the reader that 
this somewhat counterintuitive outcome is not a prob- 
lem, but rather a common feature in statistics in multi- 
dimensional spaces. 



V. RESULTS 

The f{R) functions for the HS and St models in 
Egs. fT^ and (fT5|) depend on three parameters, while 
other three parameters are needed to set the initial con- 
ditions (fTU)). By adding the GRBs intrinsic scatter aGRB, 
we end up with a seven dimensional parameter space to 
be explored. To do this efficiently, we use a Markov Chain 
Monte Carlo code which maximizes the likelihood 
along a chain with ^ 400, 000 points. The constraints on 
the parameters are then obtained by cutting out the first 
30% of the chain; we thus skip the burn-in phase and 
thin the chain to reduce spurious correlations. 

Before discussing the results, we warn the reader that, 
in the context of Bayesian statistics, the best fit model, 
i.e, the set of parameters p maximizing C{p), represents 



A. The HS model 

The HS f{R) functional expression depends on three 
dimensionless positively defined parameters. While it is 
reasonable to expect that n is not a large quantity, noth- 
ing can be said a priori on the order of magnitude of 
(ci, C2). Indeed, should both of them be very large, then 
eq. (fT^ reduces to /(i?, ci 3> 1, C2 ^ 1) ~ — m^(ci/c2) so 
that this case provide an effective cosmological constant. 
In contrast, should C2 be very small, we get f{R) ^ i?" 
thus recovering power -law corrections to the Einstein - 
Hilbert Lagrangian which are known to provide acceler- 
ated expansions for n ^ 2. We must therefore explore 
a huge range for (ci,C2); we thus skip to logarithmic 
units and get constraints on (logci,logC2) by running 
the MCMC code. 

The best fit model turns out to be : 



^ Note that a similar term is also present for SNela, but it is es- 
timated to be (Jint = 0.15 and yet included into the error ai 
provided in the Union dataset. 



nM = 0.282 , h = 0.703 , qa = -0.67 , 



n = 4.26 , logci = -0.53 , logc2 = -8.39 , 
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X 


XBF 


{x) 


^med 


68% CL 


95% CL 




0.282 


0.282 


0.282 


(0.268, 0.296) 


(0.256, 0.309) 


h 


0.703 


0.700 


0.700 


(0.692, 0.708) 


(0.685, 0.716) 


go 


-0.67 


-0.63 


-0.62 


(-0.70, -0.55) 


(-0.88, -0.45) 


n 


4.26 


2.52 


2.36 


(1.84, 3.24) 


(1.44, 4.21) 


log Ci 


-0.53 


1.15 


1.00 


(-1.07, 2.96) 


(-2.31, 7.12) 


log C2 


-8.39 


-7.64 


-7.65 


(-9.23, -6.11) 


(-9.90, -5.27) 




0.41 


0.29 


0.29 


(0.17, 0.41) 


(0.03, 0.50) 
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TABLE I: Constraints on the HS model parameters. 



-2 




■1.0 -0.9 -0.8 -0.7 -0.6 -0.5 

^0 



with a best fit GRB intrinsic scatter acRB = 0.41. The 
overall quality of the fit may be seen by looking at Fig.[T] 
and quantified by considering the following estimators : 

XsNeia/d.o.f. = 1.03 , xhRB/d.o.f. = 1.17 , u;m = 0.139 

so that we can safely conclude that the HS model is in 
very good agreement with the data. A cautionary note 
is in order here concerning the values reported above. 
The MCMC code maximizes the likelihood (|14p which 
is strictly not the same as minimizing either XsNeia '-'^ 
Xgrb- Moreover, from a statistical point of view, the sig- 
nificance level of the above reduced values can not be 
estimated from the usual tables since these standard re- 
sults do not take into account systematic errors or intrin- 
sic scatter. As such, both x'sNeia/d-o-.f- and XoRB/d-O-f- 
must be considered only as a useful tool to quantify the 
agreement with the data, but should not be overrated. 

The constraints on the single parameters are summa- 
rized in Table I where we give the mean and median 
values and the 68% and 95% confidence limits. While 
the best fit solution for {ilM,h,qo) is quite close to the 
median values, this is not the case for (n, logci, logC2), 
with the best fit solution for n lying marginally outside 
the 95% CL. However, according to the Bayesian philos- 
ophy, one must take the constraints in Table I as the final 
outcome of the likelihood analysis. It is, however, worth 
noting that setting all the parameters to their median 
values gives : 

XSNela/d.O.f. = 1.04 , xlRB/d.O.f. = 1.45 , = 0.139 

so that the quality of the fit is still not unreasonably bad. 

In order to further investigate the viability of the 
model, we can compare the constraints on {^^m, Qo) with 



FIG. 2: Likelihood countours in the plane (go,logci). Red, 
yellow and green shaded regions refer to the 68, 95, 99% con- 
fidence regions respectively. 



other results in the literature.^ However, the matter den- 
sity parameter Hm is always estimated by fitting a given 
model to a certain set of data so that a straightforward 
comparison may be biased by the different theory we are 
considering. Therefore, we simply note that the values in 
Table I are in very good agreement with typical estimates 
M, US from previous analyses of comparable datasets. 

Something more interesting can be said on the deceler- 
ation parameter qq. Indeed, the problem of model depen- 
dent estimates can now be avoided by resorting to cos- 
mographic analyses based only on the Taylor expansion 
of the scale factor. By using this approach, Cattoen and 
Visser ^STJ have found values between qq = —0.48 ± 0.17 
and qo = —0.75 ± 0.17, depending on the details of 
the method used to fit the SNLS dataset. A similar 
analysis, but using a different and smaller GRBs sam- 
ple, enabled Capozziello and Izzo [3^ to find values be- 
tween go = -0.94 ± 0.30 and qo = -0.39 ± 0.11, still 
in agreement with our estimates. A different approach 
has been instead adopted by Elgar0y and Multamaki 
j39| which advocate a model-independent parametriza- 



* We no longer consider anymore the Hubble constant h because 
this is typically marginalized over when fitting Hubble - diagram 
data because of the degeneracy with the (unknown) SN absolute 
magnitude. Such a degeneration is partially broken here thanks 
to the use of two diff'erent distance indicators and the Gaussian 
prior from the HST Key Project; however, we prefer not to dis- 
cuss the corresponding constraints because other model indepen- 
dent estimates in the literature (coming from, e.g., time delays in 
multiply lensed quasars) are affected by too large uncertainties. 
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tion of q{z). Depending on the SNela sample used and 
the parametrization adopted, then best fit values for qq 
range between —0.29 and —1.1 still in good agreement 
with our constraints in Table I. 

While the value of qo tells us that the model is nowa- 
days undergoing accelerated expansion, it is also impor- 
tant to check that this cosmic speed up ends well before 
the epoch of structure formation. To do so, for each 
point along the Markov chain, we compute the deceler- 
ation parameter q{z) and solve the equation q^zx) = 
with zt usually referred to as the transition redshift. We 
then make a histogram of the zt values and use this as 
a proxy for its probability distribution so finally estimat- 
ing: 

(zt) ^ 0.52 , ZT,med = 0.51 , 



68% CL : (0.44,0.58) , 95% CL : (0.37,0.70) . 

These constraints may be compared with previous re- 
sults. For instance, by using the Gold SNela sam- 
ple and linearly expanding q{z), Riess et al. [33 found 
Zt = 0.46 ± 0.13 in agreement with the updated result 



ZT 



0.49: 



-0.07 



obtained by Cuhna [40| which use the 



Union sample. Although there is a very good agreement, 
we nevertheless warn the reader that the estimate of zt is 
strongly model dependent. To be conservative, we may 
therefore only conclude that the transition to a decel- 
erated epoch in the HS model takes place early enough 
not to conflict with the growth of structures and galaxy 
formation. 

Although no previous analysis of the HS model has 
been performed, it is worth considering the constraints on 
the f{R) parameters (n, logci, logC2). The first striking 
result is the very high value of logci with logC2 being, 
in contrast, very low. Unless we are in the early universe 
when R/m? also takes a very high value, the bulk of 
the models with (logci,logC2) in Table I, the HS gravity 
Lagrangian may be approximated as : 



R+ f{R) ^R-m^ci ( 
V 



with the low value of m compensating for the large ci 
so that the two terms R and i?" are comparable. From 
Table I, we also get rt ^ 2.5 so that the likelihood analysis 
is selecting a subclass of the HS model that behaves as 
if the gravity Lagrangian were R ^ {R/m?)'^-^ in the late 
universe. 

In a sense, such a result could be anticipated by re- 
membering the history of fourth order gravity theories. 
Indeed, as quoted above, the first successful! f{R) model 
was proposed by Starobinsky in the '80s to give an ac- 
celerated expansion by adding a quadratic term to the 
Einstein - Hilbert Lagrangian, i.e. setting f{R) oc R^ . 
Indeed, the HS f{R) reduces to the Starobinsky infla- 
tionary model for n = 2 and C2 = with ci weighting the 
relative importance of the two terms R and R^ . Indeed, 



by looking at the likelihood contours in the (gojlogci) 
plane shown in Fig. [21 one can see that, for a given logci 
value, there are two go solutions, depending on the n be- 
ing larger or smaller than 2. For n > 2, the solution with 
the smaller go provides a better fit to the data. In this 
n regime, one gets that the more weight one gives to the 
corrective term, the more accelerated is the present day 
expansion. It is therefore not surprising that our likeli- 
hood analysis has indeed converged not too far from the 
inflationary f ^ R? model. 

We stress, however, that, even with ci ^ 1 and C2 ^ 1, 
the HS model remains fundamentally different from the 
Starobinsky inflationary model. Indeed, the accelerated 
expansion is now obtained at the present day rather than 
in the early universe. For very high z, Rjim? has be- 
come so large that the term C2(i?/m^)" in the denomi- 
nator of Eq. (|12p is dominant and we recover the limit 
/(i?) ~ — TO^(ci/c2), i.e. an effective cosmological con- 
stant. Moreover, it can be shown that, although large, 
this constant term is nevertheless negligible with respect 
to R so that we fully recover the GR Lagrangian and we 
do not thus violate the constraints from BBN and the 
abundance of primordial light elements. 



B. The St model 

As for the HS case, the f{R) expression of the St 
model still is a function with three parameters, namely 
(n, A, i?*). It is, however, convenient to reparameterize 
the model in terms of dimensionless quantities. To this 
aim, we first define : 



£ — i?*/i?o I Ae// — —\R^,/QH^ 



(22) 



with i?o the present day scalar curvature so that Eq. ([T 
becomes : 



/(i?) = -&HlKff 



(23) 



Concerning the values of the two parameters (e, Ag^j), 
we can only have some hints on what their order of 
magnitude can be. To this aim, we first note that 
R-k/R{) = R{zi,)/RQ, with an unknown reference red- 
shift. A plot of R{z)/Ro for a fiducial ACDM model 
shows that this ratio can easily take very large values 
even at intermediate redshift. Needless to say, the St 
model is not the ACDM one, but we can anticipate that, 
in order to fit the data, its expansion rate will not dif- 
fer much so that R{z)/Rq will take similar values. As a 
result, we therefore end up with e = R^^/Rq taking high 
values so that we skip to log e as a parameter to be con- 
strained by the MCMC analysis. To get an idea of the 
range for Ag//, we note that, in the early universe, the 
Lagrangian reduces to : 

R + fiR)^R- 6H^Aeff = i? - 2Aoo 
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Best fit St model superimposed to the data. 
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0.283 


0.278 0.277 (0.260, 0.295) 


(0.247, 0.308) 


h 




0.704 


0.705 0.705 (0.694, 0.716) 


(0.686, 0.728) 






-0.78 


-0.79 -0.78 (-0.97, -0.61) 


(-1.16, -0.52) 


n 




1.44 


— — < 0.79 


< 1.94 


loge 




1.56 


1.03 1.05 (0.53, 1.48) 


(0.24, 1.94) 


logAe// 


2.86 


2.68 2.67 (1.51, 3.91) 


(0.89, 4.42) 


CGRB 


0.38 


— — < 0.23 


< 0.49 



TABLE IL Constraints on tlie St model parameters. 

SO that Ae// = Aoo/3iJg. Since we want to recover the 
GR in this hmit, we just have to ask that R/hoa ^ 1- 
However, this ratio can easily be very small even if K^o 
is quite high so that we end up with Ae// varying over 
a wide range. Therefore, we skip again to logarithmic 
units taking logAg// as a model parameter. 

Running the MCMC gives us as best fit parameters : 

= 0.283 , h = 0.704 , go = -0.78 , 

n = 1.44 , loge = 1.56 , logA^// = 2.86 , 

with a best fit GRB intrinsic scatter (Jgrb = 0.38. The 
values of the reduced and of lum 

XSNela/d.O.f. - 1.03 , xhRB/d.O.f. - 1.22 , iUM = 0.140 , 

clearly show that the best fit St model is in very good 
agreement with the data as can also be appreciated by 
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FIG. 4: Likelihood contours in the (go,logAe//) plane. Red, 
yellow and green shaded regions refer to the 68, 95, 99% con- 
fidence regions respectively. 

looking at Fig. [3] Table II reports mean, median and 
confidence ranges for the individual parameters. An im- 
portant caveat is in order here for the parameters n and 
(^GRB- Indeed, the Markov chain tends to drift towards 
negative values of n which are a priori excluded in or- 
der the St model to pass the theoretical constraints (|lip . 
Therefore, the histogram of n values is clearly cut by 
this a priori assumption and, although formally mean 
and median values may be computed, they are not reli- 
able. We thus give only upper limits on n because the 
data are unable to constrain this model parameter. A 
similar problem also takes place for the intrinsic scatter 
cgrb which must be positive by definition. Nevertheless, 
the chain drifts towards the border probably because of 
a degeneration with n, with low n preferring low scatter. 
We therefore only give upper limits on ugrb too. 

The comparison between Table I and II shows that the 
values of the three parameters {0.m, h, go) are consistent 
with each other. There is actually some tension for the 
deceleration parameter go with the St model favouring a 
greater speed up, but the good overlap of the confidence 
ranges makes this difference statistically irrelevant. We 
do not discuss the comparison with previous results in 
the literature but we refer to what we have said above 
for the HS model. As a remark, however, we note that 
the agreement of the parameters (JIm, go) is not com- 
pletely unexpected. Indeed, these quantities set the ini- 
tial conditions pU|) so that this is only telling us that the 
auxiliary functions {yH,yR) are very close to each other 
for z — 0. Since both the HS and St models are designed 
to fit the same present day data, it is not surprising that 
they are quite similar today and lead to the same values 
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of the parameters which set the initial conditions. 

It is, on the contrary, somewhat surprising that, 
notwithstanding the initial larger acceleration, the tran- 
sition redshift (estimated with the same procedure used 
above) is lower than in the HS case, being : 



0.44 



0.43 



68% CL : (0.37,0.52) , 95% CL : (0.33,0.63) . 

Although formally in marginal disagreement at the 68% 
level with the results for the HS model, the two sets of 
constraints are however not too different. In particular, 
the St zt value remains in good agreement with the pre- 
vious estimates in the literature. 

Finally, even if we cannot put strong constraints on 
(n, loge, log Ae//), we briefly discuss their values. First, 
we note that, since the term /{eRqY quickly starts 
dominating over unity, at intermediate z, we get : 



f[R) ^ ^QH^Aeff 



R 

eRo 



1 



Let us suppose, for a moment, to drop the assumption 
n > 0. Indeed, for n < 0, the above f{R) is approximated 
by the sum of a (negative) cosmological constant and a 
quadratic R^ correction modulated by the value of loge. 
Once again, we therefore recover a Lagrangian similar 
to the inflationary Lagrangian of Starobinsky, since the 
negative cosmological constant soon becomes subdomi- 
nant. This qualitative discussion helps us understanding 
why the chain drifts towards negative n which have been 
excluded, in order to give Ea. p3)) the correct hmit for 
R ^ oo. For n > 0, the code is forced to look for an 
accelerating solution in this regime. Since, again, the 
term R^ / [eRqY quickly overcomes the unity, the above 
approximation of the St f{R) still holds; we may easily 
note that the dominant correction is of the form 1 /i?" so 
that the best fit n = 1.44 makes the model similar to the 
1 / R proposal firstly introduced as a fourth - order grav- 
ity motivated alternative to scalar field dark energy. As 
a final remark, we look at Fig. 2] which shows that go is 
correlated with the weighting parameter logAe// which 
here plays the same role as logci for the HS model. 



C. The effective dark energy EoS 

The impact of f{R) on the expansion history can be 
alternatively discussed by resorting to the effective dark 
energy equation of state (hereafter, EoS). Indeed, from 
the point of view of the background evolution, f{R) mod- 
els are equivalent to a cosmological scenario made out of 
dust matter and an effective dark energy term with an 
EoS given by : 

"2 dhiE{z) VLM{l + zf^ 



1+Weff{z) = 



3 rfln(l + z) 

^Mil + zf 



E^{z) 



1 - 



E^{z) 



(24) 



0.5 
0.25 


-0.25 
-0.5 
-0.75 
-1 



FIG. 5; Effective EoS for the HS (red solid) and St (blue 
dashed) f{R) theories setting the model parameters to their 
best fit values. 



so that the dark energy density parameter reads : 



1 - n 



M 



E^iz) 



■ exp 



1+Weff{z') 



1 



dz' 



(25) 



Fig. [5] shows the effective EoS for both the HS and St 
models setting the parameters to their best fit values. 
It is impressive to see how closely the two EoS are with 
Wef f {z) being almost exactly the same for the two models 
for z > 2. Note that here we have not plotted the full 
redshift range up to the last scattering z ~ 1100 to make 
easier to see the details at low z. However, when z ^ 0, 
we have checked that Weff{z) converges to Weff — 0, 
i.e. the effective dark energy behaves as a matter term 
in the early universe and scales with the redshift in the 
same way. This is still an expected consequence of the 
constraints (fTTj) ensuring that f{R)/R for z ^ oo. 
In such a limit, we therefore recover a GR universe with 
only matter. Hence, we get E'^{z ^ 1) ~ riM(l + z)^ 
which, inserted into Eq. ([M|) leads to Weff{z ^ 1) = 
whatever is the f{R) model considered. 

The HS and St effective EoS are different in the low z 
regime. Indeed, for the best fit models, we find : 

We ff{z = 0) = -1.09 , dweff/dz{z = 0) = 0.15 , 

for the HS model and : 

Wef f{z ^0) = -1.19 , dweff/dz{z = 0)^0.69 , 

for the St case. For both models, the present day val- 
ues of Weff and its redshift derivative disagree with 
the ACDM expected values, Weff{z = 0) = — 1 and 
dwef f I dz{z = 0) = 0. Actually, we must consider 
the Bayesian constraints obtained from evaluating these 
quantities along the chain. For Weff{z — 0), we find : 

{Weff{z = 0)) = -1.0b , [Weff{z ^0)]„,ed = -1.04, , 



68% CL : (-1.10,-0.98) , 95% CL : (-1.30,-0.89) 
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for the HS model, while for the St we get : 

K//(2 = 0)) = -1.19 , [We//(z = 0)],„ed--1.18 , 

68% CL : (-1.35,-1.03) , 95% CL : (-1.50,-0.95) . 
Defining wi ~ dwef f /dz{z — 0), we get : 

(Wi) = 0.21 , Wi^rned = 0.15 , 

68% CL : (-0.01,0.35) , 95% CL : (-0.16,1.14) , 
for the HS model, while in the St case it is : 

(Wi) = 0.83 , Wi.,ned = 0.69 , 

68% CL : (0.22,1.20) , 95% CL : (-0.03,4.05) . 

As it is apparent, an effective cosmological constant is 
reasonably consistent with these constraints. The present 
day value of the EoS is indeed close to —1, while its 
first derivative, although being not null, is quite small 
and consitent with zero within the 95% confidence range. 
However, we stress that neither the linear fit w{z) = 
Wjj+wiz nor the widely used Chevallier - Polarski - Linder 
|43| ansatz w{z) = wq + Waz/{1 + z) provide a good 
approximation over the redshift range probed. 

As a final remark, let us look again at Fig.O For both 
f{R) theories, the effective EoS takes a positive value 
during the epoch of galaxy formation {z 2) and over 
a wide redshift range in the matter dominated epoch. 
However, the effective dark energy density floE is very 
small during this period so that the deceleration param- 
eter stays almost constant to q{z) ~ 1 for z > 3 as for 
a matter only universe. It is, therefore, likely that the 
departures from Weff = have essentially no impact on 
the background when struture formation takes place. 

D. Local tests of gravity 

The f{R) functional law for both the HS and St mod- 
els has been formulated in such a way to pass the con- 
straints (fTTj) . However, this does not mean that there 
are no detectable deviations from the GR metric on So- 
lar System and galactic scales whatever the model pa- 
rameters are. Indeed, the HS and St models correctly 
represent two different classes of f{R) theories contain- 
ing, as particular cases, models which are able to evade 
the local tests of gravity. It is therefore interesting to ex- 
plore whether the cosmologically selected models belong 
to this subclass. A key role will be played by the value of 
l/flol = \df /dR\R=Rg. By evaluating this quantity along 
the chain, we get : 

(log = 3.87 , (log |/fl,o|)™ed = 3.81 , 
68% CL : (0.95,6.34) , 95% CL : (-0.17,10.01) 



for the HS model and 

(log l/flol) =0.16 , (log |/™|)med = -0.03 , 

68% CL : (-0.20,0.62) , 95% CL : (-0.32,1.50) 

for the St model. 

In order to see how these values compare with the con- 
straints from the local tests of gravity [4l| , we follow [23| 
who have explicitly considered the case of the HS model. 
They have demonstrated that, in order to be consistent 
with the constraints on the PPN parameter 7, one has 
to impose the following condition : 



/flo < 74(1.23 X 10^)"-! 



0.13 



(26) 



Such a constraint thus depends on the value of the HS 
model parameters and could, in principle, be included to 
limit the volume of the parameter space to be explored. 
However, we have preferred to let the MCMC code be free 
of searching the full parameter space since (|26p is actually 
a quite weak prior. Indeed, for the best fit model, this 
reduces to log|/j^o| < 13.4, while the best fit value is 
log Ifnol = 5.2 so that the test is easily passed. 

Unfortunately, we are unable to impose a similar con- 
straint to the St model. To this aim, one should first 
solve the field equations for a static spherically symmet- 
ric source matching the inner and outer solutions and 
taking care of a model for the Sun mass density profile. 
In [20], it is argued that, for every f{R) model, the PPN 
parameter 7 may always be expressed as,: 



7 - 



1 = -- 



2M, 



Meff+Mtot 

with Mtot the total solar mass and 



M, 



eff 



1^ 



R/n'^y'^dr 



with i?(r) and p{r) the scalar curvature of the local met- 
ric as function of the radial coordinate r and p{r) the 
mass density profile. In order to pass the Solar System 
constraints, M^ff should be as low as possible which can 
be accomplished by making R{r) follow K'^p{r) so that a 
chamaleon effect takes place. Since the St model has been 
designed in such a way to develop such an effect, we are 
confident that the Solar System constraints are fulfilled 
even if we cannot make any quantitative estimate. 

Finally, let us consider the galactic scales where the 
following constraint 



|/i?o| < 2 X 10- 



300 km/s 



(27) 



where Vmax is the maximum rotation velocity of the stars, 
has been advocated in [S^] as a f{R) indepedent condi- 
tion to be fulfilled in order to have no deviations of the 
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gravitational potential from the Ke pler ian 1/r profile on 
the galaxy scales. We remind that [23| look for a model 
describing the accelerated expansion of the Universe, but 
still require a substantial dark matter content. However, 
the leading term i?" of the Lagrangian of our best fit 
models induce a power law deviation from the Newtonian 
gravitational potential which was shown to describe the 
rotation curves of dwarf galaxies without the need of dark 
matter [i^l- Therefore, it is not surprising that our cos- 
mologically motivated constraints on log |/flol are defini- 
tively larger than the upper limit (log |/flo I) maa; — ~6 
for a Milky Way like galaxy. Although the result of the 
Lagrangian used by [42| cannot be directly extrapolated 
to the HS and St models, it is however suggestive that 
the violation of Eq. ([27|) is an expected outcome of the 
region in the parameter space preferred by the data. 

Rotation curves of spiral galaxies are an excellent 
probe of the gravitational potential so that we can rely 
on them to judge whether any deviation from the clas- 
sical Newtonian result is detected. Indeed, the observed 
flatness of Vc{r) in the outer regions, well far away from 
the edge of the visible disc, clearly indicates that some- 
thing unusual is taking place. If one assumes a priori that 
the potential remains Keplerian even on these scales, the 
only solution is to invoke the presence of a halo made 
out of dark matter particles (which have not yet directly 
observed). In contrast, one can reject the dark halo hy- 
pothesis and thus interpret the flateness problem as an 
evidence for deviations from the Newtonian potential. 
Should this interpretation be correct, one can try to fit 
the rotation curves with a modified potential coming out 
of the low energy limit of a f{R) theory, which has in- 
deed been done with success [4^]. As such, the fact that 
both the cosmologically selected subclasses of the HS and 
St models violate the constraint (^7)1 should not be con- 
sidered as a serious problem. Actually, giving off the 
dark halo in favour of a modified potential could lead to 
a different problem. In such a case, indeed, the matter 
density parameter should reduce to the baryons only one, 
i.e. Qm = — 0.04 in disagreement with the values in 
Tables I and II. A possible way out could be to invoke 
massive neutrinos [45| or a dark matter component only 
present on cluster scales. 

As a final remark, however, it is also possible that the 
constraint ([?7|) is overrestrictive or not valid at all. In- 
deed, Eq. ([7f|) has been obtained assuming a spherical 
galaxy described by an NFW [4J] profile. Such a profile 
is motivated by numerical simulations which assumes the 
validity of GR and hence a Newtonian gravitational po- 
tential. Therefore, testing the validity of GR based on a 
model which yet assumes this is the case is somewhat a 
contradiction. Moreover, the derivation in [20'| assumes 
that the galaxy is an isolated system, while galaxies ac- 
tually reside in groups or clusters. Should a long range 
fifth force be present, its effect must be taken into account 
when deriving a constraint in Ea. (j27p . As a consequence, 
an analytical computation which takes care of all these 
details is likely to be impossible and one should resort to 



N-body simulations in a f{R) cosmological background, 
which are yet unavailable. Because of these problems, we 
warn the reader against considering the violation of (j27p 
by the HS and St f{R) models as a definitive evidence 
that these theories should be ruled out. 



VI. CONCLUSIONS 

As soon as the observational and conceptual problems 
related to the cosmological constant and other dark en- 
ergy scenarios became pressing, modification of the grav- 
ity sector of Einstein field equations immediately ap- 
peared as an interesting alternative explanation of the 
observed cosmic speed up. Fourth -order gravity theo- 
ries then appeared as the most immediate generalizaton 
of Einstein GR since they just encoded all the deviations 
into a single analytic function f{R). As the problem of 
acceleration was solved in this framework, a new problem 
came out, namely how to choose a functional expression 
for f{R) which is not only able to speed up the expansion 
(there were too many actually!), but also not to violate 
the local tests of gravity and turn off its effect in the early 
universe where GR indeed correctly appears to work. 

From a mathematical point of view, one has to look for 
a f{R) expression satisfying the constraints ((TT|) which 
is indeed what the HS and St models efhcicntly do, pos- 
tulating that f{R) is given by Eqs. p^ and p^ . respec- 
tively. Our aim here was then to test whether these two 
well behaved models actually fit the data which suggest 
the accelerated expansion of the Universe. To this end, 
we have considered the background evolution as probed 
by the Hubble diagram by using both SNela and GRBs 
to cover the redshift range (0.01,6.6). As a first encour- 
aging result, it turns out that both the HS and St models 
are in very good agreement with the data. Moreover, the 
predicted values of the matter density parameter , the 
Hubble constant h, the deceleration parameter go and the 
transition redshift zt nicely compare with previous esti- 
mates in the literature. As a further positive outcome, 
the cosmologically selected subclasses of the general HS 
and St parameterizations easily pass the constraints on 
the 7 PPN parameter and reduce to GR in the early 
universe, as can also be seen by noting that the corre- 
sponding effective EoS and density parameter vanish at 
the redshift of the last scattering surface. Some unpleas- 
ant tension with the constraints on the deviation from 
the Newtonian potential on galaxy scales is possible, but 
the way these constraints are derived put serious doubts 
on their validity hence decreasing the significance of the 
problem. 

While this manuscript was near completion, a paper 
was posted on the arXiv [i^ where the authors try 
to constrain the HS model parameters using the Union 
SNela sample, the Hubble expansion and age data [47| 
and the Baryonic Acoustic Oscillations [1^. However, 
their results may not be compared straightforwardly to 
ours in Table I. Notwithstanding the different dataset 
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used, the main difference relies in how they reparam- 
eterize the HS model. While we have directly used 
(n, logci, logC2) as parameters, they prefer to use the 
set {^M ,n, ffio) with the effective matter density 
parameter and (ci,C2) estimated by: 

ci „ 1 - i^M 
— ~ 6 — = , 

C2 riM 



4 ri yciM J 

Finally, they have limited a priori their exploration of the 
parameter space to models with |//jo| < 0.1, while we 
have shown here that the best fit is obtained for much 
larger values. In principle, we can solve the above rela- 
tions to get fluj for our best fit model. But this is actually 
not possible since such relations have been obtained as- 
suming that Ifml <^ 1 which is not true in our case. Such 
a condition was also imposed in order to have We/ / {z) not 
departing too much from a cosmological constant, but we 
have shown here that this is an unnecessary assumption. 
Indeed, for both the HS and St models Weff{z) is close to 
the A EoS over a very limited range, but nevertheless we 
are able to fit the data since the impact of the effective 
dark energy term becomes quickly subdominant during 
the matter dominated era no matter what its EoS is. 

As Fig. [5] shows, although the f{R) functional expres- 
sion is different, both the HS and St best fit models have 
a very similar effective EoS and hence the predicted dis- 
tance modulus is almost the same. Moreover, the values 
in Tables I and II tell us that the SNela and GRBs Hub- 
ble diagram is unable to put strong constraints on the 
model parameters. 

It is worth wondering how the situation can be im- 
proved. As a first attempt, one can try to extend the 
redshift range probed by relying on the use of the so- 
called distance priors 0, which provide a useful set of 
distance-related quantities that summarize the informa- 
tion contained in the CMBR anisotropy spectrum. Al- 
though widely used as observational constraints in the 
recent literature [l^l, it is worth stressing that they are 
derived by assuming a fiducial ACDM model. It is then 
argued that the mean and covariance matrix of the dis- 
tance priors parameters do not change when the model 
space is enlarged, i.e. the choice of the fiducial model 
does not impact the estimate of the priors. Moreover, 



one also assumes that the posterior probability from the 
CMBR anisotropies and galaxy power spectra is correctly 
described by the distance priors, i.e. no information is 
lost when dropping the full dataset in favour of the sim- 
plified one. While both these hypotheses have been ver- 
ified for dark energy models |5l| . such an exercise has 
never been performed for f{R) theories. Therefore, to 
be conservative, we have preferred not to use them in the 
present analysis. One can, however, argue that including 
the distance priors can narrow down the constraints on 
some of the model parameters, but it is likely that the 
degeneracy between the HS and St models is not broken. 
Indeed, the distance priors probe the universe mainly at 
the last scattering redshift zls — 1100 thus complement- 
ing SNela and GRBs which cover the range (0, 7). How- 
ever, in the early universe, both the HS and St models 
converge to GR so that they are likely to be too similar 
to be discriminated by probing this redshift regime. 

Major improvements in the constraints and in the pos- 
sibility to discriminate not only between the HS and St 
models, but, more generally, among dark energy and 
f{R) theories are expected from the analysis of the 
growth of perturbations. Indeed, even if one can tailor 
the f{R) parameters in order to closely mimic the same 
expansion history of a given dark energy model, the evo- 
lution of the density perturbations can be rather different 
(27I [5^ . [ssj . In particular, this has a stron g im pact on 
both the power spectrum and halo statistics [5J] and the 
weak lensing signals [s^: it is thus possible to compare 
predictions with data and severely constrain the viability 
of f{R) theories. 

It is this combination of extended Hubble diagrams 
(made out of SNela and GRBs) and structure growth 
probes (such as galaxies power spectrum and cosmic 
shear) that will finally tell us whether the observed cos- 
mic speed up has been the first evidence of a new fluid, as 
mysterious as fascinating, or of new physics in the gravity 
sector, as unexpected as challenging. 
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